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1 Introduction 



The lattice regularization of QCD is a powerful tool to non-perturbatively 
study QCD in the low energy region. The numerical implementation in- 
volves however a finite lattice spacing a which has to be removed in the 
continuum limit. The rate of approaching the continuum limit will depend 
on the details of the lattice formulation. A systematic way to reduce the 
discretization effects order by order in a is the Symanzik improvement pro- 
gramme [1, 2, 3] for on-shell quantities [4, 5]. In the case of Wilson fermions, 
Sheikoleslami and Wohlert [6] have shown that for reducing the lattice arte- 
facts from 0(a) to 0(a 2 ) only one additional dimension five operator in the 
Lagrangian is needed. To achieve this acceleration of the continuum limit 
non-perturbatively the coefficient c sw of the corresponding operator has to be 
determined in numerical simulations. In the quenched case [7], the ALPHA 
collaboration has found [8, 9, 10] that the non-perturbative result for c sw 
deviates significantly from the one-loop perturbative value [6, 11]. Further- 
more, the effect of two species of dynamical fermions on c sw was also studied 
by the ALPHA collaboration and the difference to the quenched case was 
clearly visible [12]. The effect of a third flavor was studied by the CP-PACS 
and JLQCD collaborations [13] with the result that c sw (<7q) is not very much 
affected by it and their result is very close to the two flavor values of the 
ALPHA collaboration which can be seen in Figure 4. Our aim in this pa- 
per is to calculate c sw non-perturbatively for four flavors in the Schrodinger 
functional scheme. 

The paper is organized as follows. First we want to give a brief reminder 
of 0(a) improved lattice QCD and briefly discuss the improvement condition 
for c sw . After discussing some features of the algorithmic implementation we 
will give a summary about the simulation parameters and the raw data. The 
procedure how to determine c sw from the raw data will be described next. 
Finally we will give the conclusions. 

2 Improvement condition 

Our starting point for 0{a) improved lattice QCD is the fermion action 

S f = a 4 ^(x)[£ + m ]#r), (1) 

X 

where a is the lattice spacing and mo is the bare quark mass. The matrix D 
is the Wilson-Dirac operator [14] 

D = \ [( V* + V M ) 7m - aV* V M ] (2) 
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with the lattice covariant forward and backward derivatives V,, and V*. The 
leading order lattice artefacts in on-shell quantities which are calculated with 
this action are linear in a. However, the leading order discretization effects 
may be canceled by adding the so-called Sheikoleslami Wohlert term to the 
action [6] 

-^improved -D ~T~ C sw (J^ U J~ ^y. yo) 

The lattice field strength tensor is defined as in [15] and = f [7^, 7^]. 
The coefficient c sw in eq. (3) is a function of the bare coupling g and if it is 
chosen properly, Amproved becomes the on-shell 0(a) improved lattice Wilson- 
Dirac operator. For a complete cancellation of the 0(a) effects' 1 in correlation 
functions, the local composite fields that enter also have to be improved [15]. 
In our case, such composite fields are the isovector axial current A a ^{x) and 
the pseudo-scalar density P a (x). It turns out that only the isovector axial 
current needs an improvement because there is no dimension four operator 
with the same behavior as P a under the symmetries of the lattice theory. 
The 0(a) improvement requires the combination 

(Ai);^A; + a-c A ^ + d,)P a (4) 
where A a ^ and P a are given by 

Afa) = ^(xh^syV^) (5) 
P a (x) = ^(x) 7 5y^(x), (6) 

9/j,, <9* are the forward and backward difference operators and r a are Pauli 
matrices acting on one pair among the four degenerate flavors. The im- 
provement coefficient c A is known in perturbation theory [16] and from non- 
perturbative determinations for Nf — [7] and Nf — 2 [17, 18]. As explained 
in [7, 12], we introduce the unrenormalized PCAC quark mass 

( v _ |(jo + do)f A {x ) + c A adtd f P {x ) 
m[X0) ~ 2/ P (xo) [V 

where the correlation functions f A and fp contain and P a and are given 
by (2.1) and (2.2) of [7]. A second mass m' can be defined in the same way as 
eq. (7) but with the primed correlation functions f' A and fp ((2.5) and (2.6) 



d This refers to the massless theory which we consider in connection with the 
Schrodinger functional renormalization scheme. 
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of [7]). The unprimed and primed correlation functions are related to each 
other by a time reflection in the Schrodinger functional. Since the boundary 
conditions 

U{x,k)\ xo=0 = exp{aC fc }; C k = — diag(-vr, 0, tt) (8) 

bL 

i 

U(x,k)\ X0=T = exp{aC[}; C' k = — diag(-57r, 2vr, 3tt) (9) 

bL 

are such that Ck and C' k are not the same, fx and fx also differ. Since with 
PCAC we have, however, inserted an 'operator identity', all m(xo),m'(yo) 
differ only at the level of 0(a 2 ) effects in the improved theory. We could 
hence for some choice, such as Xq = yo = ^, impose m — m' = as one condi- 
tion for the proper choice of c sw and c\. Because the coefficient ca is a priori 
not known it is advantageous however to first eliminate this parameter and 
define a quark mass M independent of ca which agrees with the quark mass 
m up to 0(a 2 ) effects. For this purpose we name the partial contributions 

= (gj + w»[-) m 

4/p(a?o) 
adod fp(x ) 
= 2/ P (,o) (U) 

and rewrite the mass m as 

m(x ) = r(x ) + c A s(x ). (12) 

With an analogous definition for m', a quark mass M can than be written in 
the following way 

w/ x r \ t \ m (yo) - m '(yo) ( \ i \ r (yo) - r '(yo) 

M{x , yo) = m(x ) - s(x ) — — = r(x ) - s(x )^^ 

(13) 

In this combination, which in the improved theory differs from m by 0(a 2 ) 
only, ca drops out. Now, we define M'(x ,2/o) analogously and could require 
that the difference 

AM (^-T, ]t) = M (\t, -t) - M' (It, \t) (14) 



4 '4 y \4 '4 / V 4 4 

has to vanishes for our value of c sw . The choice (xo,yo) = (|T, \T) is one 
possible choice [7] for the argument of AM. For the quark mass M itself we 
choose (xo,yo) — {\T, \T) [7]. In order to reproduce the tree level value of 
c sw exactly for finite a, we finally impose the improvement condition 

AM = AM (0) (15) 
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where AM' ' is the tree level value of perturbation theory in the O(a) im- 
proved theory. For L/a = 8 one finds for example [7] 

aAM (0) |m=o, Csw =i = 0.000277. (16) 

As discussed in [19] , improvement coefficients possess a unique perturbative 
expansion but non-perturbatively they are themselves ambiguous by cutoff 
terms and thus depend on the choice of the improvement conditions. Of 
course, for QCD these uncertainties amount to cutoff effects beyond the order 
that is improved, 0(a 2 ) in the case at hand. In principle one then has to 
determine improvement coefficients for one fixed set of conditions as functions 
of go at constant physics, i. e. fixing all scale ratios except a that shrinks 
with go . For the Schrodinger functional this would in particular require 
constant L/tq as go is lowered. For practical reasons we fix however a/L 
instead and refer the reader to sect. (1.2.4.1) of [19] for a detailed discussion. 
The replacement of zero by the small tree-level value on the right hand side 
of eq. (16) guarantees that our definition has the correct limit for g — > 0. 



3 Simulations 

Our simulations for Nf — 4 are based on an adaptation of TAO codes — suit- 
able for APE computers [20] — used earlier by the ALPHA collaboration for 
Nf = 2 studies. An ordinary HMC algorithm [21] has been implemented with 
symmetric even-odd preconditioning [22, 23, 24] and the Sexton- Weingarten 
integration scheme [25]. Mass preconditioning [26, 27] was not enforced as 
we expect that the gain for the Schrodinger functional with the parameters 
envisaged here would not be so significant [28] . For the sake of convenience, 
from here on, we will set the lattice spacing a to one. 



3.1 Algorithm 

The fermionic determinant for an even number of flavors can be represented 
by pseudo-fermion fields fa, fa,i = 1, . . . , Nf/2. The partition function then 
reads for Nf — 4 

Z = J DUD(j){DfaD(j)\D(j) 2 e~ 5G - E -i.2^ (QQt)_1 ^ (17) 

where Sq is Wilson's plaquette gauge action [14] and Q is related to the 
Dirac matrix M by 



Q = 75 M = 75 



M ee M eo 
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The key idea of the even-odd preconditioning [22, 23] is to divide the lat- 
tice sites into even and odd sites according to the sum over the coordinates 
Yj%=o x p- ^ ft i s even for the lattice site x, the site is called even otherwise 
it is called odd. Following this strategy, the Dirac matrix M decomposes into 
the block form in eq. (18). The components are given as follows [29] 

M X;X > = (1 + T XiX )S x>x i 

- « ^ [(1 - 7m)^)<W,*' + C 1 + TaOOJO* - P)5x-H,x>] • (19) 

The matrices T r r 



-c^na^F^x) (20) 



T = 

vanish if the improvement coefficient c sw is set to zero and the submatrices 
M ee and M QO then become equal to the unit matrix. In our case however, 
we consider the 0(a) improved Sheikoleslami-Wohlert action where the co- 
efficient c sw is determined non-perturbatively for Nf = 4. Furthermore, the 
submatrices of M possess the following properties 

Ml = M ee Ml = M QO 

Ml = 75 M e 7 5 Ml = 7 5 M eo75 . 

The origin of the algorithmic acceleration by even-odd preconditioning lies in 
the factorization of the determinant. In this context, two possibilities appear. 
We may either factorize out only M ee or M OQ (asymmetric), or second, we 
extract both factors (symmetric). The derivation of both versions is very 
similar. The key point is the calculation of a determinant of matrices like 
eq. (18) using its Schur complement. For a general matrix which consist of 
submatrices, the determinant can be calculated in the following way 



det ( C D ) = det ^ AD ~ ACA ' lB }- 



(22) 



In our case, the determinant of Q = 75 M leads to 

det Q = det{M ee } det{M GO } det{<5} (23) 

for the symmetric even-odd preconditioning, where 

Q = 75(1 - MoeMgg 1 M eo ) • (24) 

In our numerical implementation, we only use the above described symmetric 
even-odd preconditioning but there is no fundamental problem to implement 
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the asymmetric even-odd preconditioned version. However, it should be kept 
in mind that the authors of [30] found that the performance of HMC with 
symmetric even-odd preconditioning is roughly 30% better than the HMC 
algorithm with asymmetric even-odd preconditioning. Due to the decompo- 
sition in eq. (24), only the odd components (f> of the fields (f> appear. With 
symmetric preconditioning the partition function for four flavors now reads 

Z = J DU D(j) ol D(f>l 2 D(j) o2 e - SG[u] - s ^ [u] ~ s ^ [u ^'' t> ^ (25) 

with the gauge part Sg[C/] an d 

S dct = N { [lndet{M cc } + lndet{M 00 }] (26) 

3.2 Simulation parameters and raw results 

The simulations were performed in the Schrodinger functional scheme [8, 9, 
10] with periodic boundary conditions for the spatial extension and Dirichlet 
boundary conditions in the temporal direction. This means the phase 9 was 
here set to zero in all runs. The data were obtained on hypercubic Euclidean 
16 x8 3 lattices. The 0(a) improvement of the Schrodinger functional requires 
additional improvement terms at the boundaries. However, since PCAC is 
an operator relation, these terms are irrelevant for a correct determination of 
c sw . Nevertheless, we have chosen them as follows. The pure gauge part of 
the action acquires a weight w(p) = c t (g ) of time-like plaquettes p attached 
to the boundary planes for which we have inserted the 2-loop value [31] 

ct^o) = 1 + [-0.08900 + 0.0191410JV f ]^ 

+ [-0.0294 + 0.002iV f + 0N?]g% , (28) 

and for the fermionic improvement coefficient Ct(<7o) [16] we take the 1-loop 
perturbative value [32] 

Ct(0 O ) = l-O.O1795sg. (29) 

In all our simulations the trajectory length was kept fixed to one. We chose 
the intervals in f3 similar to [12]. The mean acceptance rate was around 90% 
and we performed 4400 trajectories per value of j3 and c sw on average. A large 
part of our computations ran on APEmille machines with 128 processors 
each. For some values of (5 and c sw we also used apeNEXT crates with 256 
processors. The summary table of the measurements can be found in the 
appendix. 
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4 Determination of c, 



The numerical procedure for the determination of c sw involves the following 
main steps 

1. Compute AM and M for several k at fixed g\ (or (3) and c sw . Then 
interpolate linearly in M to find AM at vanishing quark mass M — 0. 

2. For fixed #q, repeat step 1 for several values of c sw and find c* w which 
solves eq. (15) by a linear fit in c sw . 

3. Repeat the preceding steps for a sufficient range of g% and fit these 
data with an appropriate function to represent the smooth functional 
dependence of c sw on g%. 

This procedure is computer time demanding because for each value of (3 and 
for each value of c sw of Table 2, we would need several runs, at least three, 
for interpolating AM in M to M = 0. To save computer time, we modified 
this method slightly. As discussed in the determination of c sw for Nf = 0, 2 
[7, 12] , the weak dependence of AM on M holds also at Nf = 4. For one set 
of parameters (3 and c sw , we checked the dependence explicitly (Fig. 1). Since 

x 1 o" 3 



f I 5 

-(5.1 -0.05 M 0.05 0.1 0.15 

Figure 1: Mass dependence of AM at (3 = 5.0 and c sw = 2.4. The errors of 
M are smaller than the symbol sizes. 

AM depends weakly on M, we contented ourselves with determining AM 
for some \M\ < 0.03 [12] and used these values of AM as an approximation 
for AM at M = 0. A typical result is shown in Figure 2. We performed 
simulations from f3 — 12 to /3 = 5.0. For each /3, we calculated the observables 
AM and M at least for three different values of c sw in such a way that the 
condition |M| < 0.03 held and that AM had a change of sign. The linear 




S 




Figure 2: Determination of c sw at (3 — 5.4. The desired value c* w is located 
at the point where the fit curve (solid) passes through zero (dashed). 



interpolation AM = s (c sw -c* w ) + AM' ' yields the desired values c* w shown 
in Table 1. We also tried to go below (3 = 5.0 but at (3 = 4.8, we were not 



p 


c* 


P 


c* 


12 


1.1429(39) 


6.0 


1.463(19) 


9.6 


1.1895(62) 


5.7 


1.554(17) 


7.4 


1.2955(76) 


5.4 


1.583(25) 


6.8 


1.3375(94) 


5.2 


1.614(28) 


6.3 


1.389(12) 


5.0 


1.717(31) 



Table 1: Results of the linear interpolation 



able to locate a significant sign change of AM in our data and therefore after 
some attempts, we decided to stop searching. 

The CP-PACS and JLQCD collaborations computed c sw for Nf = 3 in the 
Schrodinger functional setup of lattice QCD with the plaquette gauge action 
[13]. They found that the result for three flavors is very close to the two flavor 
result [12]. In addition they calculated c sw with four flavors for (3 = 9.6 and 
found c* w = 1.1954(48) in good agreement with our c* w = 1.1895(62) . 

After obtaining the proper values c* w which satisfy the improvement con- 
dition eq. (15), we want to represent and interpolate our data by a simple 
Pade formula, appropriate for the achieved precision, which also incorporates 
the known 1-loop perturbative result. The solution that we want to advocate 
here for Nf = 4 is 

c ( , 2) _ l-0-1372, 2 - 0.1641^ + 0.1679^ 

CsA9o)~ l- 0.4031s 2 <s <1.2. (30) 
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This curve appears as Fitl in Figure 3. The two other lines Fit2 and Fit3 




9l 

Figure 3: Comparison of different Pade- approximation formula for our data. 

include one and two more powers in the numerator. The deviation of the 
different fit formulae in the range g% e [0, 1.2] (/? = [12.0,5.0]) is negligible 
and beyond g$ = 1.2, Fitl and Fit2 are almost the same down to j3 = 4.5 but 
Fit3 deviates slightly. Our non-perturbatively determined formula eq. (30) 
for c sw with four flavors is valid down to (3 — 5.0 (g$ = 1.2) but may perhaps 
be used to (3 = 4.5 within a small uncertainty. 



5 1.5 
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N f =0 




N =2 




N f =3 




- N =4 


/ '/ X 


pertur. 




^^^^^ 



0.2 0.4 0.6 0.8 1 1.2 1.4 

9l 



Figure 4: Summary plot of all known c sw ((/q, iVf) for the plaquette gauge 
action. 

10 



To conclude in Figure 4, we juxtapose our new data and fit formula at 
N { = 4 to those known for N { = 0, 2, 3. 

5 Conclusions 

In the present paper, we have performed simulations for the calculation of 
the improvement coefficient c sw for four flavors of Wilson fermions in a range 
of (3 > 5.0 and give eq. (30) as a suitable parameterization of our data. We 
compared three simple-minded different Pade-approximation formulae for 
our data and could show that a small extrapolation of c sw beyond (3 = 5.0 to 
(3 = 4.5 with the main formula eq. (30) may still be acceptable. An immediate 
use of c sw will be the determination of the non-perturbative running of the 
Schrodinger functional coupling for Nf = 4 massless flavors, extending the 
present knowledge which comprises N{ = [7], iVf = 2 [12] and Nf — 3 [13]. 
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Appendix 

Table of data: 





K 


Csw 


M 


AM 


12 


0.130280 
0.129897 
0.129449 


1.028654 
1.128654 
1.228654 


0.000224(75) 
-0.001510(74) 
-0.001371(82) 


0.001735(89) 
0.000468(84) 
-0.000827(86) 


9.6 


0.131516 
0.131164 
0.131164 


1.140488 
1.170488 
1.250488 


-0.00004(11) 
0.00595(11) 
-0.00853(11) 


0.00077(12) 
0.00053(11) 
-0.00040(12) 


7.4 


0.134626 
0.133753 
0.132989 
0.132349 


1.163222 
1.263222 
1.363222 
1.443222 


-0.00196(20) 
0.00072(20) 
0.00018(20) 
0.00027(20) 


0.00169(19) 
0.00069(15) 
-0.00036(16) 
-0.00162(17) 


6.8 


0.135638 
0.135082 
0.134896 
0.133813 
0.133056 


1.209613 
1.299613 
1.309613 
1.409613 
1.509613 


0.00128(30) 
-0.00642(32) 
-0.00333(29) 

0.00320(29) 
-0.00141(30) 


0.00169(26) 
0.00033(26) 
0.00079(24) 
-0.00000(24) 
-0.00204(22) 


6.3 


0.137098 
0.136018 
0.135028 
0.134028 


1.239058 
1.339058 
1.439058 
1.550580 


-0.00129(38) 
0.00142(36) 
0.00140(36) 

-0.00251(31) 


0.00199(24) 
0.00077(27) 
-0.00075(27) 
-0.00115(24) 


6.0 


0.138358 
0.136669 
0.134375 


1.250000 
1.387912 
1.587912 


-0.00090(85) 
0.00214(47) 
0.00677(43) 


0.00270(30) 
0.00022(48) 
-0.00123(34) 


5.7 


0.140327 
0.138229 
0.137008 
0.135685 
0.133940 


1.250000 
1.387912 
1.487912 
1.587912 
1.754350 


-0.0082(16) 
0.00034(63) 
0.00081(55) 
0.00564(63) 
-0.00004(38) 


0.00369(38) 
0.00107(43) 
0.00076(48) 
-0.00032(37) 
-0.00156(28) 


5.4 


0.141417 
0.139111 
0.137815 
0.135028 
0.133775 


1.307912 
1.487912 
1.587912 
1.787912 
1.907912 


-0.00032(77) 
-0.00712(89) 
-0.00762(70) 
-0.00117(69) 
-0.00854(39) 


0.00258(36) 
0.00112(55) 
-0.00056(57) 
-0.00164(47) 
-0.00202(31) 


5.2 


0.143363 
0.140628 
0.139206 
0.138147 
0.136248 
0.134556 


1.307912 
1.487912 
1.587912 
1.655891 
1.787912 
1.907912 


-0.0004(14) 
-0.00032(87) 
-0.00326(69) 
0.00162(94) 
0.00030(91) 
0.00372(62) 


0.00196(48) 
0.00076(42) 
0.00102(39) 
0.00028(47) 
-0.00084(65) 
-0.00189(38) 
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D 1 460^6 


1 307Q1 9 


n nnM (9^ 


D 009^9(69) 




D 149^4 


1 507919 


0091 (1 3) 


001 95(41 1 

\J .\J\J -LsLiOy L ±± 1 


5.0 


n i 38141 


1 78791 9 


—0 0053(1 1 


001 09(68") 


136527 


1 885463 


0009(111 


00004(55) 




0.135039 


2.000000 


-0.00826(90) 


-0.00070(59) 




0.129603 


2.400000 


0.00033(42) 


-0.00482(27) 




0.145928 


1.500000 


0.0083(54) 


-0.0002(12) 




0.142201 


1.700000 


-0.0022(31) 


-0.0000(11) 




0.138295 


1.910000 


0.0017(13) 


0.00060(62) 


4.8 


0.137971 


1.930000 


0.0031(12) 


-0.00054(48) 




0.136844 


2.000000 


0.0008(11) 


-0.00044(50) 




0.135327 


2.100000 


-0.00492(87) 


-0.00116(44) 




0.132358 


2.300000 


-0.00550(72) 


-0.00244(41) 



Table 2: Summary table of the measurements 
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